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Abstract This paper deals with the vibration of granular materials due 
to cyclic external excitation. It highlights the effect of the acceleration on 
the settlement speed and proves the existence of a relationship between set- 
tlement and loss of contacts in partially confined granular materials under 
vibration. The numerical simulations are carried out using the Molecular Dy- 
namics method, where the discrete elements consist of polygonal grains. The 
data analyses are conducted based on multivariate autoregressive models to 
describe the settlement and permanent contacts number with respect to the 
number of loading cycles. 
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1 Introduction 

Ballast materials of railway platforms exhibit complex behavior under re- 
peated loading. With the increase of vehicle speeds as well as comfort and 
safety requirements, understanding the dynamics of these materials is be- 
coming a crucial issue. UnHke highly agitated granular materials which can 
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be described with the kinetic theory, ballasted layers are generally subjected 
to dense flow where the trajectories of grains are correlated and the coUisions 
cannot be considered as randomly distributed in terms of positions and veloc- 
ities. Therefore, such a statistical approach may be inappropriate to predict 
the response of railway platforms. 

In this research work, the Molecular Dynamics method is used to sim- 
ulate granular samples made of polygonal grains under vibration. Since its 
introduction by Cundall and Strack ||9], this discrete elements method has 
proved its viability in describing several mechanisms such as granular mate- 
rials transport \2u\, mixing [3], segregation compaction [ia] etc. Coupled 
with the experimental approaches, the MD is now recognized as a fundamen- 
tal tool to inve stig ate the behavior of dense granular materials. Recently, Lu 
and McDowell [19| developed an approach based on the MD to investigate the 
permanent displacements in granular beds, made of irregular shaped grains, 
under a single loading cycle. The complex geometry was produced using 
spherical grains assemblies. Although efficient in terms of grain shapes, this 
technique neglects the inertial effects of the overlapping regions. The discrete 
element method was also adopted by Lobo-Guerrero and Vallejo [3| on rail- 
way platforms to investigate the effect of grain degradation on granular bed 
response under cyclic loading. The model took into account the rupture of 
grains using a criterion based on the loading modes and force intensities. Sim- 
ilar approaches such as the molecular dynamics method was used to simulate 
the vibration of confined granular material made of polygonal components 
[1; [23|. In this paper, the target is to relate the permanent settlement in 
partially confined samples to the loss of contacts using statistical analysis. In 
the second section, a brief description of the simulation method is suggested. 
The settlement mechanism which represents the residual displacement un- 
der sleeper is described through field cases simulation. In the third section, 
a causality analysis based on multivariate autoregressive models describing 
the settlement and contact loss is conducted. 



2 Simulation method 

The sample is composed of polygonal grains which are described with a set of 
ordered vertices (sa,i)jG[|i,6|] ■ The positions of the vertices can be evaluated 
using the distances ri and r2 as well as the orientation 6, which are shown 
in figure[l]-(a). These parameters follow uniform bounded distributions: ri ~ 

^r,„i„,r„ax]i '^2 C^[0.25ri,0.75ri], and 9 ^ C/[0:27r] ■ 

The interaction between grains takes place when an overlap is detected 
(Figure [TJ-b) . The dichotomy method is used to calculate the shortest dis- 
tance between the vertices. This leads to the edges which are candidates for 
interaction. Among this Hst of edges, the contact segment which relates the 
couple of intersection points is obtained. It defines the tangential component. 
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t, of the contact frame and its perpendicular represents the normal compo- 
nent, n. The reference point of the frame, o, is defined by the middle of the 
above mentioned intersection points. Once the contact frame, the overlap, 
and the velocities are known, the contact forces acting on a grain a can be 
written as follows: 

* \sign{F^)^iFn - ^tVt if ||Ft|| > /iF„ 

The interaction forces acting on the particles contain elastic terms; by as- 
suming uniform pressure and shear stresses along the contact surface, it can 
be shown that these elastic forces are Hnear: F^ = knUx and F^ = ktUy, 
where /c„ = ^ (i-iy'^) ^^'^ ~ et(i-iy^) • ^ ^^'^ ^ materials properties 

(Young modulus and Poisson's ratio) whereas the constants c„ and Ct can 
be obtained experimentally using a compression test. The interactions also 
enclose viscous terms denoted by the phenomenological constants 7„,t which 
can be rewritten as dissipation fraction with respect to the critical damping: 
7n,t = c(n,t\/mkn,t- Finally, the Coulomb friction is included through the 
threshold Ft < /ii^„, where fi is the friction coefficient. At the same time 
the moments are calculated using the contact frame and the contact forces. 
The discrete element simulations are carried out with the following param- 
eters: density of grains, pp = 2710kg/m^, Young modulus, E = 46.9GPa, 
Poisson's ratio, v — 0.25, friction coefHcient, /i = 0.8, viscous coefficients, 
ttn = 0.8 and at = 0.1, and grain dimensions rmin = Tmax = 5 mm (unless 
differently specified). It is worthwhile noticing that the interaction with the 
wall is described the same way as the interactions between grains. 

Using the contact forces between two grains a and (3 at the reference 
point o, it is possible to calculate the moments around the centroid Cq, using 
the identity M^.^ = Mo + i^c x oCq,, where F^. is the contact force applied on 
the grain a as expressed beforehand in equation |(T]). The procedure is used 
for each component of the sample. Once the moments and interactions are 
known, the equations of motion can be integrated with respect to time using 
a finite difference scheme according to the Molecular Dynamics 

method. 



2.1 Sample preparation, loading, and boundary conditions 

At the beginning of the process, the grains are subjected to the gravity field 
until reaching the full equilibrium. The displacements and rotations are cal- 
culated using the predictor-corrector algorithm [l|, where the contact forces, 
the moments and the body forces are taken into consideration. Once the 
equilibrium is reached, the sample - of width R = 75 mm, height H = 150 
mm and number of grains 175 - is subjected to a sinusoidal load of the form 
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F ^ Fo + AF sin ust, where uj is the circular frequency, Fq is the initial force 
which is taken equal to 0.5kN, and AF is the force amplitude. The frequency 
and force amplitude vary in such a way that the sleeper covers a wide range 
of accelerations around g (gravity). The excitation is applied on the inner 
half (0 < r < 0.5i?) and upper end {yt=o = H) of the sample through a rigid 
sleeper, unless differently specified. The sides r — R and y = represent 
the wall (container), they react to the grains actions as described before- 
hand. However, at r = 0, a symmetry condition is simulated by omitting the 
friction effect (Figure [J]- (a)). Henceforth, this loading case will be termed as 
partially confined configuration. The fully confined configuration corresponds 
to the case where the excitation is applied on the whole upper end of the 
sample. 

Under repeated loading granular materials undergo large displacements 
towards free regions. It is interesting to notice that observations of displace- 
ment fields and strains showed that granular materials deform because of 
rearrangements of the packing, rather than contact elasticity QUI. Unlike 
fully confined samples, where the settlement is mainly due to grain rear- 
rangement and overwork, the mobility of grains in partially confined granular 
materials is predominant. An experimental work performed by the authors 
14 1 on irregular ballasted samples showed that the acceleration of the 



sleeper plays a key role on the mobility of grains. In both confinement cases, 
it has been shown that the settlement increases with the acceleration. In 
addition, in case of partially confined samples, a sharp increase in terms of 
settlement speed has been noticed when the acceleration exceeds the gravity. 
At high level of agitation, the loss of contacts is frequent especially at the 
critical plane relating the edge of the sleeper to the wall as can be seen in 
figure {lb). 



2.2 Settlement in terms of acceleration and loss of contacts 

The simulations conducted herein are performed using polygonal grains in 
order to produce a shape which is similar to the micro-ballast used in the 



above mentioned study |l3l : Il4| . The obtained numerical results consist of 
axial displacement of the sleeper with respect to time. The applied forces 
are varied from 0.5 to 3 kN and the frequency is varied from 20 to 40 Hz, 
in such a way that the sleeper covers a wide range of accelerations around 
the gravity. The grains radii are distributed uniformly from rmin = 5 mm to 
rmax = 10 mm. In accordance with the experimental procedure, the settle- 
ment under sleeper, hmax is described with a logarithmic law of the form: 
hmax{N) = A + B\n{N), with respect to the number of cycles TV. This log- 
arithmic law seems to be valid for different granular materials under cyclic 
loading 0; [ill; EH| The parameters A and B depend on different phys- 
ical factors such as the degree of confinement, the applied force, and the 
frequency. These factors are independent and they can affect individually 
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the settlement speed. The experimental result revealed that the acceleration 
which depends simultaneously on the applied force and frequency is the best 
explicative quantity in terms of correlation. Moreover, it has been shown that 
the degree of confinement results in different behaviors in terms of settlement 
versus acceleration. The simulations conducted herein produce most of the 
experimental features. It can be seen that the partially confined configura- 
tions exhibit a critical transition in terms of settlement speed at around lAg, 
as can be seen in figure This transition has been observed experimen- 
tally it is probably due to the loss of contacts which can not be detected 
easily with the considered experimental setup. Numerical simulation can pro- 
vide some more information regarding the history of the granular material 
texture. 

Figure ^ shows the variation of the sleeper displacement as well as the 
number of contacts with respect to the number of loading cycles, at differ- 
ent frequencies. In these particular cases, the number of cycles is limited to 
around 25, for clarity, and the applied force amplitude is equivalent to 1.5 kN. 
It can be noticed that the settlement is much higher when the acceleration is 
beyond the gravity level (/ = 40 Hz). On the other hand, it can be seen that 
the number of contacts follows the exciting force, in terms of oscillations. 
This means that while cyclically loaded, the granular material exhibits local 
alternating opening and closure of contacts independently of the acceleration 
level. However, it can be seen that there is a difference of about 20% in terms 
of permanent loss of contacts when the sleeper undergoes high acceleration. 

Under cyclic loading of partially confined sample, the grains fiow towards 
less loaded regions as can be seen in figure ^h) . More interestingly, it can be 
seen in figure ^a) that during the settlement process, contact openings occur 
at a specific region between the front of the sleeper and the wall frontier. In 
the same region there is a loss of density as can be seen in figure ^h) . 

This critical and localized behavior can be observed during the settlement 
process independently of the number of cycles. In the following section, it will 
be shown that there is a causality relationship between the loss of contacts 
and the settlement speed. 



3 Causality relationship between contacts loss and settlement 
speed 

In order to find out a logical relationship between the responses of gran- 
ular materials under vibration, statistical analysis are necessary since the 
solution is not analytically determinist. In this section, the objective is to 
investigate the causality between the settlement velocity and the loss of con- 



tacts. Therefore, we adopt Granger's approach which was introduced in 
econometric in order to forecast possible relationship between discrete time 
series. However, this approach can be applied for physical systems providing 
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accessible responses with respect to time [5|. The model questions whether 
the prediction of a given variable is improved by taking into account its own 
history and the history of another variable. The original model concerns only 
stationary time signals (time series where the first and second moments are 
independent of time) or Hnear (time series with a stationary rate of change) . 
More recently, it has been shown that it is possible to conduct causal- 
ity testing on non stationary time series if (i) they can be approximated by 
tendency functions and (ii) the estimation errors are stationary. 



3.1 Causality measurement 

The physical quantities of interest are the settlement under sleeper, h(t), and 
the number of permanent contacts, z(t). In order to illustrate the causality 
analysis, a granular sample consisting of polygonal grains of diorite with an 
average size of 5 mm is considered. The applied force in this particular case 
is of frequency 60 Hz and amplitude IkN. The time series are extracted 
from the numerical results by averaging over the loading cycles as follows: 

nT 

z{n) z{t) >~ ^ / z{t)dt and h{n) —< h{t) >, where T is the loading 

(n-l)T 

period (Figure (H). The dynamic relationship between the above mentioned 
variables are then described with a multivariate autoregressive model of the 
order / as follows: 



h{i) = ^'{i) + akh{i -k) + J2 bkz{i - k) + e,,(i) 

T T (2) 

z{i) = z'^{i) + Ckh{i — k) + dkz{i — fc) + ez(i) 

k=l fc=l 

where h^'{i) = a + /3Ln{i) and z'^{i) = x + 7* ^ire the tendency functions 
with respect to the number of cycles, and e are the error terms, which are 
assumed to be independent. 

In accordance with the Granger's method, a comparaison in terms of 
accuracy is conducted between the above description ^ and the following 
model: 



h{i) — h^{i) + ^ akh{i - k) + eh{i) 

k=l 

I 

z{i) = z''{i) + ^ dkz{i - k) + ez{i) 



T (3) 



k=l 



where h^{i) = a + (3Ln{i) and z^{i) = x + ji are tendency functions with 
respect to the number of cycles, and e are the error terms which are assumed 
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to be independent. The main idea is to compare the accuracies of the two 
models. If the description ^ improves the estimation of the physical quantity 
h as compared to the description that means that z causes h, since z 
is an expHcative variable of h. The non causality hypothesis is expressed by 
Hq : bi = b2 = ... = bi = 0. Henceforth, the term "i/o" will be used to refer 
to the above mentioned hypothesis. When this hypothesis is vaHd, the first 
term of the system ([2]) is reduced to the first term of ([3]). Therefore, it is 
possible to calculate the causality measure of z and h based on the above 
mentioned autoregressive models as follows: 



where N denotes the length of the times series, and represent the 
auto- variances of Cs and es where s — h,z. Under the nullit y hy pothesis Hq, 
the measure J-z->h has an asymptotic x^(2^) distribution In order to 
reject the hypothesis Hq with a preselected risk level a, the calculated 
should be higher than the value of the standard distribution of parameters 
a and I). 

3.2 Estimation coefficients 

The causality measure between the considered signals can be calculated using 
the error terms and eg, where s = h, z. Therefore, it is necessary to estimate 
first of all the coefiicients of the suggested models and deduce the error terms. 
This task can be accomplished using the recursive least square method 
In this section, we adapt our description to the algorithm by expressing the 
signals under consideration as follows: 



where d{i) is an output corresponding to h, z or both of them as will be seen 
later on, A is a vector which encloses the unknown estimation parameters 
and u{i) represents the estimation error at the step i. In the case of model 
(121), these vectors can be written as an input a;*^*' [1, Ln{i), hii — 1), .., h{i — 
z(i — 1), .., z{i — I)] , an unknown A = [a, /?, ai, .., a;, x, 7, 6i, .., bi] and 
error terms u{i) — th{i)- Similarly, in the case of model ([31), these vectors can 
be written as a;'*^ = [1, Ln{i), h{i — 1), .., h{i — I)], A — [a, P, oi, .., a/] and 
u{i) = eh{i). In order to estimate the unknown parameters of the model |(5]), 
the least squares criterion is used: 




(4) 



d{i) = A'x'-'^ + uii) 



(5) 



n 




(6) 



Solving the problem consists in minimizing the quantity TZ with respect to 
the unknown vector A. This leads to a recursive algorithm (appendix), which 
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consists in starting from an initial set of parameters A ~ 0, an initial matrix 
Q = / of size I X I and a real constant << A < 1. At each iteration n, the 
calculation steps read: 

. _ \-\Qx 

u{n) = d{n) ~ x^A 

A = A + ku{n) ^ ' 

Q ^ {Q - kx*Q) 

where the vector k is termed as the Kalman gain. Knowing the parameters 
of the models ^ and |[3]), it is possible to plot out the estimations of the 
settlement and permanent contacts number as shown in figure 

Once the error terms are obtained, the causality measure can be calcu- 
lated using equation In the particular case considered in this study, the 
discrete time signals h and z are of length N = 1000 cycles and the order of 
autoregression is of ^ = N/100. The risk level of rejecting the hypothesis Hq 
is a = 1%. At this risk level, and for a number of parameters I, the calculated 
theoretical value of X^(20 is 37.52, however, the measure of causality equals 
590.85. Therefore, it can be concluded that the number of permanent con- 
tacts loss is a significant expHcation variable of the settlement under sleeper. 
This approach can be applied for different loading cases. Table H]) shows that 
the causality direction remains valid for different frequencies and amplitudes 
of loading. 



4 Conclusion 

In this study, numerical simulations of a partially confined granular material 
under vibration are presented. It has been shown that the settlement process 
is characterized by a fiow of materials towards less loaded regions. It has also 
been noticed that there is privileged regions of contacts loss. Furthermore, it 
has been proved that the loss of contacts causes the settlement, for different 
loading cases. The residual displacements which take place under dynamic 
loading at different frequencies and amplitudes depend on several factors 
such as material properties, grain shape, degree of confinement etc. In this 
study, we selected the diorite and the polygonal shape because this type of 
material is widely used in railway platforms. Moreover, we concentrated on 
the loss of contacts as an explicative variable of settlement. As a perspective 
of the suggested analysis, it would be interesting to investigate all the physical 
factors which may infiuence the settlement speed, such as the acceleration 
and elastic defiection. It is also possible to extend the suggested procedure for 
more complex grain shapes. For instance, as long as the settlement depends 
on the mobility of the material, non convex grains may have an important 
effect on the settlement. Actually, unHke convex grains (circular, polygonal, 
elliptic etc.) where the contact are binary, non convex grains can be connected 
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with more than a single contact. Assuming equivalent sizes of contact areas, 
the dissipation should be higher in case of non convex grain. In addition, the 
contact openings are more difficult to take place, therefore the settlement 
should be different. This particular point will be the subject of a future 
contribution by the authors. 



A Recursive Least Squares 

Minimizing the least square criterion leads to the following equation: 

n n 
2 = i = 

Let p'"^ = ^^^pA"-'x(i)xW, iiM = X^J^_oA"-'[xWa;«']. Using the Woodbury 
identity of matrices (Kima et Bennighof [13]) one can obtain the following relation- 
ship: 

qM ^ - A-ife("'a;("'*Q("-i' (9) 

where Q*-"' = ii~ andfe'"' = — r-n — — r^-^-r . Introducing this expression 

in the equation A'"' — Q'"'p*"\ leads to a solution of the problem in a recurrent 
formula: 

/i'"' =h("-i' +fe("'it(n) (10) 
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1.5 
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Table 1 Causality measure for different excitations 
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Fig. 1 (a) Typical grain characterized by ri, r2, and 9, (b) Contact detecting 
using the dichotomy method. 
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Fig. 2 (a) Boundary conditions and applied force (b) Residual displacement 
of grain (the colorbar represents the intensity of residual displacement in the y- 
direction in cm). 
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Fig. 3 Settlement speed with respect to acceleration. 
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Fig. 4 Response of the sample in terms of settlement {hmax) and number of 
contacts {Zr>At) with respect to the number of cycles at different frequencies. 
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Fig. 6 Variation of the settlement and of the number of permanent contacts with 
respect to the number of cycles. 
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Fig. 7 Application of the models ([2|) and ([3]) on the signals of settlement and 
permanent contacts. The subscripts "in", "1", and "2" denote the MD method cal- 
culations, the estimation using the first model and the estimation using the second 
model, respectively. 



